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ABSTRACT 

We introduce a new method for stacking voids and deriving their profile that greatly 
increases the potential of voids as a tool for precision cosmology. Given that voids are 
distinctly non-spherical and have most of their mass at their edge, voids are better 
described relative to their boundary rather than relative to their centre, as in the 
conventional spherical stacking approach. The boundary profile is obtained by com¬ 
puting the distance of each volume element from the void boundary. Voids can then 
be stacked and their profiles computed as a function of this boundary distance. This 
approach enhances the weak lensing signal of voids, both shear and convergence, by 
a factor of two when compared to the spherical stacking method. It also results in 
steeper void density profiles that are characterised by a very slow rise inside the void 
and a pronounced density ridge at the void boundary. The resulting boundary density 
profile is self-similar when rescaled by the thickness of the density ridge, implying that 
the average rescaled profile is independent of void size. The boundary velocity profile 
is characterized by outflows in the inner regions whose amplitude scales with void size, 
and by a strong inflow into the filaments and walls delimiting the void. This new pic¬ 
ture enables a straightforward discrimination between collapsing and expanding voids 
both for individual objects as well as for stacked samples. 

Key words: cosmology: theory - dark matter - large-scale structure of Universe - 
methods: data analysis 


1 INTRODUCTION 

Cosmic voids represent a potentially powerful tool for mea¬ 
suring the cosmological parameters and probing the nature 
of dark energy (e.g. Li 2011; Bos et al. 2012; Lavaux & Wan- 
delt 2012; Sutter et al. 2012b; Cai et al. 2014b, a; Hamaus 
et al. 2014a). Most cosmological constraints are derived from 
the structure and dynamics of voids, which are a probe of 
modified gravity models (Li et al. 2012; Clampitt et al. 2013; 
Cai et al. 2015; Barreira et al. 2015) as well as of the na¬ 
ture of dark matter (Hellwing et al. 2010; Yang et al. 2015; 
Massara et al. 2015). In the former case, void profiles are sen¬ 
sitive to the presence of a fifth force, which, while screened 
in higher density regions, can be large in voids. Such a force 
leads to emptier and larger voids, due to the faster evacua¬ 
tion of matter from low density regions (Peebles & Nusser 
2010; Clampitt et al. 2013). In the latter case, replacing 
cold dark matter by warm dark matter or including mas¬ 
sive neutrinos would lead to less evolved voids, and hence to 
shallower density profiles. 

Up to now, the density and velocity structure of voids 
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has been studied through the use of spherical profiles mo¬ 
tivated by the fact that stacking many voids results into 
spherically symmetric structures (e.g. van de Weygaert & 
van Kampen 1993; Padilla et al. 2005; Ceccarelli et al. 
2013; Ricciardelli et al. 2013; Hamaus et al. 2014b; Na- 
dathur et al. 2015). But individual voids are distinctly non- 
spherical. While the simple picture of an expanding under¬ 
density in a uniform background suggests that voids should 
become more spherical as they evolve (Icke 1984), in real¬ 
ity, voids are not isolated and this simplified picture does 
not hold. There are two major factors that affect the evo¬ 
lution of voids. Firstly, contrary to the case of collapsed 
structures, void evolution is strongly affected by the tidal 
held of the surrounding distribution of matter (Platen et al. 
2008; van de Weygaert & Platen 2011). Secondly, as voids ex¬ 
pand, they are squeezed by neighbouring voids. These effects 
lead to present-day voids that have highly complex shapes 
(Platen et al. 2007, 2008; Neyrinck 2008; Sutter et al. 2012a; 
Nadathur & Hotchkiss 2014). 

The diversity of void shapes makes the traditional stack¬ 
ing procedure suboptimal for extracting cosmological infor¬ 
mation. Simply put, the cosmological constraints are de¬ 
rived by comparing the density inside voids with that at 
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their boundaries. For example, in some modified theories 
of gravity the inner regions of voids are emptier than in the 
standard cosmological model, with the evacuated matter de¬ 
posited at the void boundaries. Stacking randomly oriented 
voids of various shapes leads to an overlap of the voids in¬ 
ner regions and boundaries. This “blurring” decreases the 
density contrast between the inner and outer parts of voids, 
leading to a lower signal. In addition, there is ambiguity in 
the definition of the void centre used for spherical stacking, 
with different choices resulting in different density profiles 
(e.g. Nadathur & Hotchkiss 2015b). 

In this work, we introduce a new method of both mea¬ 
suring void profiles and stacking voids by taking into account 
their shape. In contrast to the spherical method, we propose 
that void profiles should be measured with respect to the 
void boundary. This leads to a much sharper distinction be¬ 
tween the inside, boundary and outside of voids, resulting 
in at least two major gains. Firstly, it leads to a better un¬ 
derstanding of the structure and dynamics of cosmic voids 
enabling a closer comparison with analytical theories of void 
evolution. Secondly, it increases the stacked lensing signal 
of voids, which is the best probe for measuring void density 
profiles (Higuchi et al. 2013; Krause et al. 2013). 

The paper is organized as follows. In Sec. 2 we outline 
the new method by applying it to a simplified void model; in 
Sec. 3 we describe the cosmological simulation to which we 
apply the method as well as the void catalogues we construct 
from it; in Secs. 4-6 we present the density, velocity and weak 
lensing profiles obtained using the new boundary stacking 
approach. We conclude with a short discussion and summary 
in Sec. 7. 


2 THE BOUNDARY PROFILE OF VOIDS 

Here we give an overview of the proposed method for com¬ 
puting boundary void profiles, which we illustrate using a 
simplified model of a void. We construct a void by randomly 
selecting a shape for it from a cosmological N-body simu¬ 
lation. A cross section through the boundary of the void is 
shown in the top panel of Fig. 1. For simplicity, the inner re¬ 
gion of the void is assigned constant density, l + <5i ns ide = 0.1, 
where, 5 = = — 1, denotes the density contrast. The void is 
embedded within a uniform background, 1 + <5 ou t a ide = 1, 
and the mass evacuated from within the void is deposited 
uniformly on the boundary, which is shown as a solid curve. 

Finding the spherically averaged profile involves finding 
a void center, typically the volume-weighted barycentre,and 
growing concentric shells around it. This process is schemat¬ 
ically illustrated in the centre panel of Fig. 1, where, for 
clarity, we only show a few radial shells, but, in practice, we 
employ many more such shells. The spherical profile is given 
by the mean density of matter inside each shell. Applying 
this method to our model void provides the spherical density 
profile shown in the top panel of Fig. 2. For small radial dis¬ 
tances, which correspond to shells fully enclosed by the void, 
we recover the input density value, 1 + <5 = 0.1. At larger 
radii, r ^ 9 h -1 Mpc, the shells intersect the void boundary 
giving rise to a “noisy” profile. Due to the irregular shape of 
the void, different radial shells have varying degrees of over- 




Figure 1. Illustration of the new method for measuring void 
profiles. The top panel shows the void boundary, with the actual 
void shape selected randomly from voids found in an N-body 
simulation. For simplicity, the void is assigned a constant density, 
1 + ^inside = 0.1, inside its boundaries and is embedded in a 
uniform background with 1 + <5 ou t s ide = 1- The mass evacuated 
from inside the void is deposited at the void boundary, which has 
1 + ^boundary = 30. The center panel shows the spherical shells 
around the barycentre of the void that are used for computing the 
spherical profile. The bottom panel shows lines of equal distance 
from the void boundary (thick black curve) that are used for 
computing the boundary profile proposed in this paper. 

lap with the void boundary, giving rise to “noisy” features 1 . 
These persist for as long as the shells intersect the bound¬ 
ary, corresponding to r ^ 26 h _1 Mpc, while for even larger 
radii we recover the background density. This simple exam¬ 
ple illustrates that the spherical density profile is a complex 
convolution of the shape of the void and its actual density 
distribution. 

To calculate the void profile with respect to the bound¬ 
ary of the void we compute the boundary distance, V , that 
corresponds to the minimum distance from each point to the 
void boundary (see Eq. (3) for a formal definition). The out¬ 
come is illustrated in the bottom panel of Fig. 1, where each 
thin contour line corresponds to points that are at equal 
distance from the boundary of the model void. Now we can 
calculate the density profile as a function of V by computing 
the mean density inside each shell of constant V (in prac¬ 
tice, we use many more shells than those shown in Fig. 1). 

1 In contrast to our simplified model, in real voids the mass is 
not distributed uniformly along the void boundary, resulting in 
even larger “noisy” features. 
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Figure 2. Void density profile. The top panel shows the spherical 
profile of the simple void model illustrated in Fig. 1. The vertical 
grey line marks the effective radius of the void, -R e fj, defined in 
Eq. (1). The bottom panel shows the density profile of the void as 
a function of the distance, T>, from the void boundary. For clarity, 
we define T> as having negative values inside the void and positive 
outside. The vertical grey line marks to the boundary of the void. 


To distinguish between points inside and outside the void we 
adopt the convention that V takes on negative values inside 
the void and positive outside, with T> = 0 at the void bound¬ 
ary. The resulting profile is plotted in the bottom panel of 
Fig. 2 and shows that we recover the actual input density 
distribution: 1 + 5 = 0.1 inside the void, a large value of 1 + 5 
at the void boundary due to the mass evacuated from inside 
the void, and 1 + 5=1 outside the void. 

The new void profile has two major advantages. Firstly, 
it is independent of the shape of the void. For example, dis¬ 
torting the boundary of the void in Fig. 1, while keeping 
the same density distribution inside and outside the void, 
would result in exactly the same density profile as a function 
of V. Secondly, on average, the mass displaced from inside 
the void is found at the void boundary, with the resulting 
density at the boundary being at least an order of magni¬ 
tude higher than inside the void (Sheth & van de Weygaert 
2004, hereafter SvdW, and Sec. 4). Thus, while the spherical 
profile for radial shells that intersect the void boundary is 
dominated by the density at the boundary and not by the 
density inside the void, our proposed profile naturally differ¬ 
entiates between the boundary, the inside and the outside 
of the void. 


3 VOID IDENTIFICATION 

We make use of the high-resolution Millennium cosmologi¬ 
cal N-body simulation (MS; Springel et al. 2005). The MS 
follows the evolution of cold dark matter (DM) using 2160 3 
particles, each of mass, m p = 8.6 x 10 s h - 1 M 0 , to resolve 
structure formation in a periodic cube 500 h _1 Mpc on a 
side. The MS assumes the WMAP-1 cosmogony (Spergel 
et al. 2003) with the following cosmological parameters: 
flm = 0.23, Qa = 0.75, h = 0.73, n s = 1 and as = 0.9. 

We identify voids using mock catalogues constructed 
from the semi-analytic galaxy formation model of Guo et al. 
(2011). For this, we select only galaxies with stellar masses, 
M* ^ 3.8 x 10 10 h -1 M@, such that the number density is 
n = 3.2 x 10 -3 h 3 Mpc~ 3 , similar to that of typical redshift 
surveys (e.g. Zehavi et al. 2011). These galaxies are used as 
input to the Delaunay Tessellation Field Estimator (DTFE; 
Schaap & van de Weygaert 2000; van de Weygaert & Schaap 
2009; Cautun & van de Weygaert 2011), which uses a Delau¬ 
nay triangulation with the galaxies at its vertices to extrapo¬ 
late a volume filling density field. The resulting density field 
is used as input to the void identification method. We also 
apply the DTFE method to the distribution of DM particles 
to obtain continuous density and velocity fields, which are 
used for computing the density, velocity and weak lensing 
profile of voids. Both the galaxy density field and the DM 
density and velocity fields are stored on a 1280 3 regular grid 
with a grid cell size of 0.39 h - 1 Mpc. 

The voids are determined using the Watershed Void 
Finder (WVF; Platen et al. 2007), which identifies voids 
as the watershed basins of the large scale density field, sim¬ 
ilar to the ZOBOV void finder (Neyrinck 2008). Compared 
to other methods, the watershed void finders have the ad¬ 
vantage of not imposing any a priori constrains on the size, 
shape and mean underdensity of the voids they identify (Col- 
berg et al. 2008). The WVF proceeds by first smoothing the 
galaxy density field with a 2 h _1 Mpc Gaussian filter, whose 
size corresponds to the typical width of the filaments and 
walls forming the void boundaries (e.g. Cautun et al. 2013, 
2014). This smoothing is applied in order to dilute any sub¬ 
structures present on the void boundaries (e.g. see Cautun 
et al. 2014), which could potentially give rise to artificial 
voids. The smoothed density field is segmented into water¬ 
shed basins using the watershed transform implemented us¬ 
ing the steepest descent method (e.g. Bieniek & Moga 2000). 
This process is equivalent to following the path of a rain 
drop along a landscape: each volume element, in our case 
the voxel of a regular grid, is connected to the neighbour 
with the lowest density (i.e. steepest descent), with the same 
process repeated for each neighbour until a minimum of the 
density field is reached. Finally, a watershed basin is com¬ 
posed of all the voxels whose path ends at the same density 
minimum. 

To overcome oversegmentation, the WVF joins the 
basins that share a boundary with a galaxy density, 5 g ^ 
— 0 . 8 , since such low values typically separate subvoids em¬ 
bedded within larger voids. This threshold is motivated by 
the model of an expanding top-hat underdensity for which 
shell crossing takes place at 5 = —0.8 (SvdW). This top-hat 
model has several shortcomings when compared to realistic 
voids, e.g. voids do not have initial top-hat profiles, their ex¬ 
pansion is restricted by their environment and observations 
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Figure 3. The abundance of galaxy voids in the MS. The shaded 
regions show the three ranges in effective void radius, R e g, for 
which we compute average density and velocity profiles. 

provide only the galaxy density, not the total matter den¬ 
sity, which makes the extent to which <5 g ^ —0.8 is a realistic 
threshold debatable. This step leads to the merging of only 
2% of the watershed basins and hence it has no noticeable 
effect on the profiles of stacked voids. 

The distribution of the WVF voids is shown in Fig. 3 
where the voids are characterised by their effective radius, 
R c g. This corresponds to the equivalent radius of a sphere 
with the same volume as the void, i.e. 

= in 

where Woid denotes the volume of the void. The figure 
shows that we identify a wide range of void sizes, from 

5 to 50 h _1 Mpc, with the abundance peaking at R e g ~ 
15 h -1 Mpc. For the rest of this work, we will calculate 
stacked profiles for voids in three intervals in void size cor¬ 
responding to Reg = 8 — 12, 18 — 22 and 30 — 35 /i -1 Mpc 
(shown as dark shaded regions in Fig. 3), which contain 656, 
643 and 100 voids, respectively. These intervals were chosen 
to probe a variety of void sizes, while at the same time hav¬ 
ing enough voids to provide reliable statistics. 

The abundance of WVF voids is similar to that obtained 
using the ZOBOV void finder when applied to DM trac¬ 
ers with the same number density (Nadathur & Hotchkiss 
2015b, Fig. 1), but it is a factor of two higher than when 
applying ZOBOV to the galaxy distribution (Nadathur & 
Hotchkiss 2015a, Fig. 2). The difference is likely due to 
the merging criteria employed by the two void finders (see 
Nadathur & Hotchkiss 2015a who studied the dependence 
of the void abundance on merging criteria). Regardless of 
these differences, the WVF voids have a similar minimum 
galaxy density to the ZOBOV voids (see Fig. 1 in Nadathur 

6 Hotchkiss 2015a). 

3.1 Spherical profiles 

The spherical profile of a void is computed as a function of 
the radial distance from the void centre, which we take as the 
volume-weighted barycentre. While there are other potential 


choices of void centre (e.g. see Nadathur & Hotchkiss 2015b) 
that result in slightly different spherical profiles, these dif¬ 
ferences are small when compared to the difference between 
the spherical and the boundary profile methods. The void 
centre is given by x vc = where the sum is over 

all the N voxels that are part of the void and Xi gives the 
position of each such voxel. The density of the void at radial 
distance, r, is then given as 

S(r) = Wk5k , (2) 

where the sum is over all the voxels found at a radial dis¬ 
tance, r ± | A r, with A7' the radial bin width, and 5k is the 
density of each voxel. The weights, Wk , give the overlapping 
volume between the voxel and the radial bin. This is cal¬ 
culated by generating 64 points regularly distributed inside 
each voxel; Wk is then the fraction of those points that are 
found inside the radial bin. This method of calculating void 
profiles is very similar to the VTFE method of Nadathur 
et al. (2015), except that we use a Delaunay instead of a 
Voronoi triangulation. Besides the density profile, we also 
compute the profile of the velocity component along the ra¬ 
dial direction, Vm. This is computed similarly to Eq. (2), but 
with the density replaced by the radial component of the 
velocity, vii = v ■ r/r. 

3.2 Boundary profiles 

To calculate the shape independent profile we need to iden¬ 
tify the boundary or border of each void, which is the density 
ridge delimiting the watershed basin corresponding to that 
void. In practice, we compute the void boundary as follows. 
We loop over all the neighbouring grid cells of each voxel 
that is part of a void. If one of the neighbours is not part 
of the same void, then the face connecting the two voxels 
is identified as the boundary of the void. To speed up the 
computation, each such face is stored as only one point cor¬ 
responding to the centre of the face. Finally, the border of 
the void is given by the union of all those points, i.e. by all 
the centres of the faces connecting voxels that are not part 
of the same void. This procedure can be easily expanded to 
identify the boundary of ZOBOV voids too. In this case, the 
voids are composed of Voronoi cells, not the cells of a reg¬ 
ular grid as for the WVF. The boundary of ZOBOV voids 
is given by the union of the faces of the Voronoi cells that 
connect two Voronoi cells that are not part of the same void. 

The next step is to compute the distance of each point 
from the void boundary, which in computer science is re¬ 
ferred to as the distance transform. This technique has been 
previously used to find the galaxies that are the farthest in¬ 
side voids (Kreckel et al. 2011) and to use voids for improv¬ 
ing photometric redshift estimation (Aragon-Calvo et al. 
2015). The minimum distance from the void boundary to 
a point with Cartesian coordinate, x, is given by 

+ min; [ |x — y;| ] for x outside the void ^ 

— mini [ |x — y i | ] for x inside the void 

where {yi} denotes the set of points that give the void 

boundary and | | denotes the magnitude of a vector. By con¬ 
vention, the boundary distance is negative for points inside 
the void and positive outside. One can further define a void 
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boundary distance field, which at each point in space is a 
vector of magnitude, \D\. given by 

© = © X ~ Yj , (4) 

\x-yA 

where j denotes the index of the point on the void boundary 
closest to x. The direction of © is perpendicular to the sur¬ 
faces of constant © (see bottom panel of Fig. 1) and always 
points outwards. The void boundary distance field is com¬ 
puted separately for each void, using a kd-tree constructed 
from the set of points that gives the boundary of the void. 

For each void, the boundary distance takes a minimum 
value, ©min, that corresponds to the point inside the void 
farthest from the boundary. We always have —©min ^ R e g 
(note that © m i n is negative), with equality only for spher¬ 
ical voids; on average, we find © m i n /-R e ff = —(0.62 ± 0.06) 
(ler standard deviation). For spherical voids, the boundary 
distance, V, is equivalent to the radial position, and, for 
this particular case, the spherical and boundary profiles are 
exactly the same. 

The void density profile as a function of V is computed 
similarly to Eq. (2): 

m = , ( 5 ) 

T,k W k 

but now the sum is over voxels found at a distance ©± | A© 
from the void boundary, with A© the width of the © bin. 
The weight, Wk, is given by the fraction of the 64 uniformly 
distributed points inside each voxel that are within the re¬ 
quired distance from the void boundary. As in Sec. 3.1, we 
define the velocity component along the direction of T> as 
U|| (©) = v ■ ©/©. 

3.3 Stacking 

The stacked profile is computed as an average over the indi¬ 
vidual profiles of voids in a narrow R e g range. For spherical 
profiles, we average as a function of the rescaled radial dis¬ 
tance, r/R e fi. For boundary profiles, we average over the 
individual voids at constant boundary distance, ©. In this 
latter case, the minimum distance, ©min, can differ even be¬ 
tween voids of equal radii. As a result, at very low © values 
only a subset of the voids contribute to the stacked bound¬ 
ary profile. This effect, which our calculation accounts for, 
becomes important only for the lowest © bins and can be 
easily spotted due to the large error bars associated with 
these points. The uncertainty for the two stacking proce¬ 
dures is given as the lcr interval in the distribution of the 
mean value obtained using 1000 bootstrap samples. 


4 VOID DENSITY PROFILE 

We now apply the boundary profile analysis to galaxy voids 
found in the MS. To demonstrate the power of this new 
method, we compare with the outcome of the conventional 
approach based on spherically averaged profiles. 

4.1 Individual voids 

In Fig. 4 we show the density profile of six random voids se¬ 
lected to span a wide range of sizes. These profiles and the 
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Figure 4. The density, 1 + 5, profile of six randomly selected 
voids that span a range of effective void sizes, -R e fj. The top panel 
shows the spherical profile as a function of the rescaled radial 
position, r/R e g. The bottom panel shows the boundary profile 
as a function of the distance, T>, from the boundary of the void. 


subsequent ones are computed using the full DM particle 
distribution and hence give the overdensity of matter. The 
spherical profile is shown as a function of the rescaled radial 
distance, r/R e g, which is used for determining self-similar 
and universal void profiles (Hamaus et al. 2014b; Nadathur 
et al. 2015). While there is a large variation between the 
different voids and between neighbouring radial bins, on av¬ 
erage individual voids are underdense for r/R c ft <0.5 and 
show no consistent features at larger distances. Two of the 
voids have S > 0 at their centres that can be explained ei¬ 
ther by the presence of a substructure inside the void (Beygu 
et al. 2013; Rieder et al. 2013) or by the void centre being 
close to the boundary (Nadathur & Hotchkiss 2015b). 

Compared to the spherical profile, the boundary profile 
is very different and has more features, indicative of the fact 
that, since voids have highly complex shapes, taking a spher¬ 
ical average erases or damps many features. In addition, the 
boundary profile shows a better qualitative agreement be¬ 
tween the various voids: underdense for © ^ —3 h -1 Mpc; a 
sharp density peak at the void boundary, © = 0; and close 
to mean background density for © ^ 5 /i _1 Mpc. The density 
peak at the void boundary is expected since voids identified 
using watershed-based methods are delimited by a density 
ridge. The height of this ridge is largely given by the mass 
contained in the most massive haloes, which explains the 
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Figure 5. The stacked density profile for voids in three ranges 
in effective radius, i? e jf. The two panels show the spherical (top) 
and the boundary (bottom) profiles. 


variation in height between different voids. Massive haloes 
can also be found outside the void, resulting in sporadic 
peaks in the density profile, but only very rarely inside the 
void - no such example is present in Fig. 4. The width of the 
density peak at the boundary is given by the typical size of 
the massive haloes as well as that of the filaments and walls 
that delimit the void (e.g. Cautun et al. 2013, 2014). 


4.2 Stacked profiles 

In Fig. 5 we present the mean density profiles of voids in 
three R e g bins chosen to probe a variety of void sizes (see 
Fig. 3). The spherical stacked profiles are underdense in the 
inner parts, with <5 slowly rising to a maximum at r ~ R e ft, 
followed by a gradual transition towards the average back¬ 
ground density (Hamaus et al. 2014b; Nadathur et al. 2015). 

The boundary profile paints a different picture of the 
structure of voids. In the inner most parts, T>< —4 h _1 Mpc, 
the density is very low, —0.9 ^ 5 ^ —0.5, and nearly con¬ 
stant, with only a very small increase in 8 with V. This is 
followed by a very steep rise of a density ridge at the bound¬ 
ary, which decreases nearly as fast at V ^ 0. At even fur¬ 
ther distances, the density gradually reaches the background 
value. 

The boundary density profile can be understood within 
the multiscale picture of the cosmic web. Void interiors are 
not fully empty, but instead are criss-crossed by tenuous fil¬ 



Figure 6. A simple model to understand the boundary profile. 
The thick black curves show the boundary of the central void 
and that of its neighbours, which are coloured according to their 
density, with dark and light grey showing high and low density. 
The highest density regions correspond to the intersection points 
of two or more void boundaries, with the density decreasing far¬ 
ther away. The thin curves shows contours of constant distance, 
T >, from the boundary of the central void, with two of those con¬ 
tours, T> = —5 and 5 h _1 Mpc, highlighted in cyan. The outer 
contours intersect the boundary of neighbouring voids and hence 
correspond to a higher mean density than the inner contours. 


aments and walls that become more densely packed as one 
approaches the massive structures that delimit the voids 
(Cautun et al. 2014). Thus, the mean density is expected 
to increase close to the void boundary, in accord with the 
results shown in the bottom panel of Fig. 5. Close to the 
void boundary, the behaviour is dominated by the promi¬ 
nent filaments and sheets that delimit the void and that are 
substantially denser than the tenuous structures found in¬ 
side the void (Cautun et al. 2014). The picture outside the 
void is complicated by the presence of neighbouring voids 
and their own dense ridges, as illustrated in Fig. 6. The 
density profile is not symmetric around T> = 0 since neigh¬ 
bouring voids can have different sizes, and hence different 
ridge thicknesses. In addition, the outer contours intersect 
the boundary of neighbouring voids. Due to clustering, the 
density varies along the void ridge, with higher density typ¬ 
ically associated with the intersection points of two or more 
void boundaries. The density profile is sensitive to this clus¬ 
tering, which would explain why the slope, |, is shallower 
outside the void than inside the void. 

Compared to the spherical profile, the average boundary 
profile shows smaller differences between voids of different 
sizes and is close to a self-similar profile. Before discussing 
these differences, we proceed by fitting the boundary profile 
with the empirical function: 


P = 


Pin (i + (^ - l) e (1 -q|D|) 

P °** (i+fc- 1 ) e "^) 


for V < 0 

1 

for V ^ 0 


( 6 ) 

where p = p(l + <S) is the matter density and p is the mean 
background density. The fit is a continuous function com¬ 
posed of two parts that describe the inner, T> < 0, and outer, 
T> ^ 0, mean density profiles, with p(T> = 0) = p max . 

The very interiors of the void are characterised by the 
density parameter, pi n , and by the slope a, the latter ac¬ 
counting for the fact that the density increases with V. The 
density ridge at V ~ 0 is well described by an exponential 
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Figure 7. The best-fitting function (Eq. 6) to the boundary 
density profile of voids. The top panel shows the mean density 
for voids with g = 18 — 22 h~ 1 Mpc (symbols with error bars) 
and the best-fit function (dashed line). The remaining panels show 
the ratio between the data and the best-fitting function for voids 
of different sizes. The fit was done using only data points with 
V < 10 h- 1 Mpc. 


function that takes a maximum value, p max , at V = 0. This 
ridge is not symmetric with respect to V = 0 and so we 
have two parameters in the exponential, tj n and tout, that 
give the thickness of the inner and outer void boundary, re¬ 
spectively. Just outside the void boundary, the density has 
yet to converge to the background value, so there is an addi¬ 
tional parameter, p out , to account for this effect. The 0^0 
part of the fitting function should include an additional com¬ 
ponent to account for the transition towards the background 
density at large V , but, for simplicity, we omit such a compo¬ 
nent. Our function is characterised by six parameters which 
is similar to other empirical fits to spherical void profiles: 
Hamaus et al. (2014b) proposed a four parameter fit that 
was latter extended by Barreira et al. (2015) to a five pa¬ 
rameter fit to give a better description of voids identified 
in the DM density field. Compared to the boundary profile, 
the spherical one smooths over many density features, so it 
is not surprising that the former requires more parameters. 

The upper panel of Fig. 7 shows that the empirical 
function of Eq. (6) describes, to very good approximation, 
the mean density profile. To better assess the fit quality, 
the lower panels of Fig. 7 show the ratio between the mea¬ 
sured profile and the best-fitting value for three void sam¬ 
ples. The fit matches the data well, except for a few points 
around T> ~ 0, which show a ~10% difference, and for the 
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Figure 8. The best-fitting parameters of Eq. (6) obtained from 
stacked void density profiles. The top panel show the thickness 
of the inner void boundary, t[ n , as a function of void radius. The 
remaining panels show the dependence of the other fit parameters: 
Pin, a, Pmax, Pout and i ou t as a function of / j n . The error bars give 
the ltr uncertainty. The dashed lines show that the best-fitting 
parameters follow simple relations with R e g (top panel) and ti n 
(remaining panels). 


D ^ —15 h _1 Mpc region of the largest voids, which shows 
a systematic deviation from the best-fit. 

Fig. 8 shows the best-fitting parameters and their lcr 
errors for voids of different size. These were computed us¬ 
ing the Markov chain Monte Carlo method implemented 
in the emcee package (Foreman-Mackey et al. 2013). The 
figure shows tin as a function of i? e ff and the remaining 
parameters as a function of tj n . The best-fitting parame¬ 
ters follow linear relations with ti n , which in turn can be 
parametrized as a quadratic function of This suggests 
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that the parametrization of Eq. (6) is overdetermined and 
that the number of free parameters is too large (similar re¬ 
lations between the fit parameters have been reported by 
Hamaus et al. 2014b). Eq. (6) can be rewritten by express¬ 
ing pi n , a , p m ax, Pout and tout as a linear function of tin (two 
parameters in each case) and, in turn, by expressing ti n as a 
quadratic function of i?. e ff (three parameters). This results 
in a 13 parameter function that fits in one step voids of var¬ 
ious sizes. We repeated the fit using these parametrizations 
and obtained similarly good fits. 

According to Fig. 8, void interiors are characterised by 
a nearly constant density, pi n , but by different values of the 
density slope, a, with larger voids having more slowly vary¬ 
ing density profiles. The height of the density ridge, p ma x, 
is largest for small voids since these are typically embedded 
in overdense regions. This is illustrated also by the p ou t/p 
density parameter that is larger than 1 for the smallest voids 
and that decreases with void size. The density ridge is asym¬ 
metric and is thinner inside the void, i.e. tin < tout (see the 
discussion of Fig. 6). 

We also find that the smallest voids have lower ti n values 
and larger t ou t values than the largest voids. The increase 
of tin and decrease of tout with void size can be a manifes¬ 
tation of the age characterising voids of different size. Just 
as low mass haloes, small voids are dynamically old, so the 
density ridge has been squeezed for a longer time. Larger 
voids, which originate from larger scale density fluctuations, 
have not had enough time to pile up mass at the ridge to 
the same extent as the small ones. 

4.3 The self-similarity of stacked profiles 

The boundary density profile of voids of different size is very 
similar, but not exactly the same (see bottom panel of Fig. 
5). Those differences are minimized, or even disappear en¬ 
tirely, when rescaling the inner profile by the thickness, tin, 
of the inner void boundary. The rescaled profiles are given 
in Fig. 9 which clearly shows that all voids, independently 
of their size, have a self-similar profile. To better highlight 
this, in the bottom panel of the figure we take the ratio with 
respect to a weighted mean density. This weighted mean was 
obtained by averaging, at fixed 'D/tin values, over voids of 
different sizes, with the contribution of each sample weighted 
by the inverse of its associated error. Small systematic differ¬ 
ences with void size are seen only for V ~ 0, which probably 
arise because small void are embedded in overdense regions 
while large voids are found in predominantly underdense re¬ 
gions. For the rest, all the density profiles lie on the same 
curve with less than ~5% scatter. 

The self-similar nature of boundary profiles suggest that 
voids of different sizes have, on average, the same dynamical 
characteristics. In contrast, the same self-similarity is not 
seen for spherical profiles (see top panel of Fig. 5). This could 
be due to the limitations of spherical profiles, among which, 
most importantly, is the mixing and inability to separate 
between the inside, boundary and outside of voids, as we 
exemplified in Sec. 2. This fits in with the results shown in 
the bottom panel of Fig. 5 where the differences between 
voids of various sizes are most pronounced in the boundary 
and outside regions of the voids. 

Self-similar profiles are obtained only after rescaling by 
the thickness of the inner void ridge, ti„. This suggests that 



-20 -15 -10 -5 0 


rescaled boundary distance T) / t in 

Figure 9. The self-similarity of voids. Top panel: the density 
profile, 1 + 5, as a function of the rescaled void boundary distance, 
T>/ti n , where t; n is the thickness of the inner void boundary as 
determined by fitting Eq. (6) to the density profile. The symbols 
correspond to voids of various effective radii, R e g. All voids have 
a self-similar density profile independent of Bottom panel: 

the ratio between the profiles and a weighted mean of the values 
in the various R e g bins showing that there is less than 5% scatter 
among voids of various sizes. 

the void interior knows about the boundary or vice-versa, 
and that the two evolve together. The former possibility 
seems ruled out by the simple picture of an expanding spher¬ 
ical underdensity in which the evolution of a shell of matter 
of radius, r, depends only on the mass contained within r 
(Fillmore & Goldreich 1984, SvdW, but see Ruiz et al. 2015). 

Spherical void profiles have also been claimed to be 
self-similar (Ricciardelli et al. 2014; Nadathur et al. 2015, 
e.g.), but there are contradictory results in the literature 
(e.g. Hamaus et al. 2014b; Sutter et al. 2014; Nadathur & 
Hotchkiss 2015b, this work). The self-similarity of spherical 
profiles seems to be dependent on several factors: the void 
finder, the population of tracers used to identify the voids 
and the tracers used to measure the void profile. This could 
be the case for boundary profiles too, though it is reassuring 
that self-similarity of boundary profiles has been found for 
both voids identified using galaxies (this work) and for voids 
identified in the DM density field (Cautun et al. 2015). 

4.4 Comparison to analytical predictions 

It is illustrative to compare with analytical predictions of 
void profiles, among which the isolated spherical underden¬ 
sity model is the most popular (Fillmore & Goldreich 1984, 
SvdW; see van de Weygaert & Platen 2011 for a more elab¬ 
orate void evolution model that includes ellipsoidal collapse 
and that accounts for the effect of the external tidal field). 
For this purpose, we select a top-hat spherical underdensity 
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Figure 10. Comparison of analytical and measured density pro¬ 
files of voids. The solid line corresponds to an uncompensated 
top-hat spherical underdensity that gives rise to a void with 
mean density, 1 + <5 = 0.3. The dotted and dashed curves give 
the spherical and boundary distance profiles of MS voids with 
R e fj = 18 — 22 /i _1 Mpc. The top-hat void shows a good qualita¬ 
tive agreement with the boundary distance profile of MS voids. 


that gives rise to a void of radius, R e a = 20 /i _1 Mpc, and 
density, 1 + 5 = —0.3, similar to the mean density of stacked 
MS voids with sizes, f? e ff = 18 — 22 /i _1 Mpc. While realis¬ 
tic voids do not have initial top-hat profiles, such a simple 
model captures most of the features of initial underdensi¬ 
ties representative of cosmological environments (see Fig. 3 
of SvdW). Fig. 10 shows the density profile of the resulting 
void as a function of the rescaled radial distance, r/R e g. 
The figure also shows the spherical and boundary profile of 
MS voids with sizes, R e a = 18 — 22 h _1 Mpc. To plot all 
three profiles on the same x-axis, we give the boundary pro¬ 
file in terms of the rescaled coordinate, (V + R b h)/R e s, with 
R c f f = 20 h -1 Mpc. 

The top-hat profile shows large differences with respect 
to the spherical profile of MS voids, but is in approximate 
agreement with the boundary profile of the same voids. In 
particular, the boundary profile matches the main prediction 
of the analytical model, the formation of a density ridge 
at the edge of the void. Thus, this simple model offers a 
qualitative description of the density profiles of voids, but 
only after accounting for the fact that real voids are non- 
spheric.al. 

Note, however, that there are significant differences be¬ 
tween the top-hat model and the boundary density profile 
of realistic voids, which are driven by many factors. Our 
goal is not to test the accuracy of the analytical model, but 
rather to show that such a model performs better than one 
would naively expect from a comparison to spherical pro¬ 
files. For example, the edge of MS voids contains more mass 
than the analytical prediction since the boundaries of realis¬ 
tic voids accrete matter also from outside the void (note the 
1 + 5 < 1 values of the boundary profile for rescaled radial 
positions larger than 1.3). Secondly, replacing the uniform 
top-hat underdensity with more realistic initial density pro¬ 
files results in a more gradual increase of the density ridge 
(SvdW), which is closer to the profile of MS voids. 

Fig. 10 is a first step towards testing one of the central 


assumptions of the SvdW void abundance model, which is 
that void formation is well described by the evolution of an 
isolated spherical underdensity. We have shown that we can 
find a top-hat model that qualitatively matches the mean 
density of stacked voids. It remains to be seen if the pa¬ 
rameters of this top-hat model are also the ones required to 
match the initial conditions of realistic voids. Furthermore, 
for the model to be realistic, the match should work not only 
for stacked samples, but also for individual voids. 


5 VOID VELOCITY PROFILE 


The velocity field of voids is another property that can be 
better understood by analysing boundary profiles. As for 
the density profile, we proceed by comparing the spherical 
and boundary velocity profiles. We focus on the peculiar 
velocity component, Wm , that gives the rate at which matter 
is evacuated in comoving coordinates through a surface of 
r = const and T) = const for the spherical and boundary 
profiles, respectively. Positive U|| values correspond to a net 
outflow of matter through the surface while negative values 
correspond to an inflow. 

For investigating void velocity profiles we use the same 
objects, both individual and stacked samples of voids, as 
we used when studying the density profiles in Sec. 4. Figs 
11 and 12 show the corresponding wm profiles for individual 
and stacked voids. For brevity, we focus our discussion on 
the stacked velocity profiles, with individual voids showing 
similar trends, albeit with large individual variations. 

The spherical velocity profiles show outflows from voids, 
which peak at ~0.6R e ff, and that are followed by regions 
with lower outflow velocities or even inflows. The nearly- 
linear increasing outflow for r< 0.5 R e a indicates that void 
interiors expand faster than the average universe showing 
a so-called super-Hubble outflow (Icke 1984; van de Wey- 
gaert & Platen 2011; Aragon-Calvo & Szalay 2013). For the 
boundary profile, the velocity, vy, increases until near the 
void boundary, T>< — 3 h -1 Mpc, and is then followed by a 
rapid switch from outflow to inflow. This behaviour at V ~ 0 
is consistent with infall onto the void boundary, which, given 
its high density, is the main local driver of dynamics. At fur¬ 
ther distances from the void boundary, the velocity slowly 
converges towards 0, as expected. 

Given the density profiles shown in Fig. 5, we can use 
linear theory to predict the U|| values (e.g. see van de Wey- 
gaert & van Kampen 1993), which are shown as dotted lines 
in Fig. 12. The linear predictions are given by 


V||,lm = — 


Hf M{< x ) 
Pm S(X) 


( 7 ) 


with H the Hubble factor, / ~ Q°; 55 the linear growth factor 
and p m the mean background density of matter. The symbol 
x stands for the radial distance, r, for spherical profiles and 
for the distance, V, for boundary profiles. The factor M{< 
x) denotes the mass contrast inside x and S(x) denotes the 
area of a surface of constant x. See Appendix A for details 
and for a short derivation of the relation. 

The linear theory prediction agrees with the data for 
the spherical profile, except for a few small systematic ef¬ 
fects: the velocity of small voids is overpredicted while that 
of large voids is underpredicted. These discrepancies, seen 
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Figure 11. The peculiar velocity profile of the six randomly 
selected voids shown in Fig. 4. It shows the velocity component, 
Di|, along the direction of r and T>, respectively. The two panels 
show the spherical (top) and boundary (bottom) profiles of those 
voids. 

also by Hamaus et al. (2014b), have been attributed to the 
effect of surrounding structures on void interiors (Ruiz et al. 
2015). In the case of boundary profiles, the linear theory is in 
agreement only for the void interior, i.e. T>< 0, and at large 
distances, T>> 10 /i -1 Mpc. Large discrepancies are present 
at the void boundary and just outside the void where the 
linear predictions can be off by up to 100 km/s. Such dif¬ 
ferences are not surprising since linear theory is valid in the 
regime |<5| 1. For spherical stacking, while the average 5 is 

not very small, it is below unity at every point. In contrast, 
the boundary stacking has very large values of 5, as high 
as 3 at the void edge, which explains why large discrepan¬ 
cies are seen only at, and just outside, the void boundary. 
Moreover, as shown in Fig. 4, individual voids have density 
values above unity for spherical profiles as well, so linear 
theory would break down in such cases too. The difference 
is that for spherical profiles the position of the 5 > 1 region 
varies from void to void, so departures from linear theory 
average out when stacking many such objects, whereas for 
the boundary profile the departures are always at the same 
position, V ~ 0. 

We find that the U|| peak is highest for spherical profiles 
and that the same peak is up to 20% lower for boundary 
profiles, even though in the latter case the velocity increases 


Figure 12. The peculiar velocity profile as a function of radial 
distance, r, for spherical stacking (top) and as a function of void 
boundary distance, D, for boundary stacking (bottom). It shows 
the velocity component, ny, along the direction of r and T>, re¬ 
spectively. To guide the eye, the data points are connected with 
solid lines. The dotted lines show the linear theory prediction for 
uij given the average density profiles of Fig. 5. 
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Figure 13. The boundary profiles of the velocity component 
along the radial vector, r, (solid with symbols) and along the 
boundary distance vector, T> (dotted). It shows that the outflow 
from voids is preferentially along the radial direction increasing 
until close to the void edge. 
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Figure 14. The velocity of the void boundary, U||. boundary) as 
a function of void size. Negative values correspond to contract¬ 
ing voids and positive values to expanding voids. The top panel 
shows this velocity for voids stacked according to their size, i? e fj. 
It shows the velocity at T> = 0 (solid curve) and the mean ve¬ 
locity over the interval \T>\ r 1 h~ 1 Mpc (dashed curve), which 
is more robust. The bottom panel shows the probability distri¬ 
bution function (PDF) of the ridge velocity for individual voids 
of various sizes. The distribution is very broad with each sample 
having both expanding and contracting voids. 


until close to the void edge. To explain this, Fig. 13 shows the 
boundary profile for the velocity component along the radial 
direction. For comparison, the dashed lines show the profile 
of the velocity component along 27, which corresponds to 
the solid lines with symbols in the bottom panel of Fig. 
12. For 27 < 0, the radial velocity is larger than the velocity 
component along 27 which shows that the outflow from voids 
is preferentially directed radially. Fig. 13 also shows that the 
radial velocity, when binned according to 27, increases until 
close to the void boundary and it is then followed by a very 
steep decrease at the edge of the void. This contrasts with 
the spherical profile of the radial velocity (see top panel of 
Fig. 12), which shows a peak at r ~ 0.6 R e s and not at the 
edge of the void, i.e. r ~ R e a- 

The boundary profile offers a natural way of discrim¬ 
inating between expanding and contracting voids. For ex¬ 
ample, expanding voids correspond to a positive uy value at 
their boundary since the boundary is moving outwards. The 
top panel of Fig. 14 shows the values of the velocity at the 


boundary, 27 = 0, and also the vy value averaged over the 
interval |27| ^ 1 /i _1 Mpc, with the latter being less prone 
to noise. The plot shows that, on average, small voids are 
contracting while large ones are expanding, with voids of 
R e ft ~ 25 /i -1 Mpc being at the transition between the two 
behaviours. Using the mean density of the large-scale region 
in which voids are embedded, previous studies have charac¬ 
terised the voids as under- or overcompensated (e.g. Cecca- 
relli et al. 2013; Hamaus et al. 2014b; Nadathur & Hotchkiss 
2015a). For example, the smallest two void stacks in Fig. 
12 are overcompensated while the larger voids are slightly 
undercompensated. This distinction can be determined, for 
example, using the sign of the radial velocity at r> 1.5f? e ff 
(see top panel of Fig. 12), with positive values corresponding 
to underdense regions and vice versa. Combining this with 
our analysis of the void boundary dynamics, we find that 
overcompensated voids are predominantly contracting while 
the undercompensated ones are predominantly expanding. 

Using the boundary profile one can determine even for 
individual voids if they are expanding or contracting, as we 
show in the bottom panel of Fig. 14. For example, while 
most small voids are contracting, there is also a significant 
fraction that are expanding. Similarly for the largest voids: 
while most are expanding, there are large contracting voids 
too. Thus, expanding and contracting voids cannot be differ¬ 
entiated using just their size, R c a, and additional void prop¬ 
erties need to be considered (see e.g. Nadathur & Hotchkiss 
2015a). 


6 WEAK LENSING FROM VOIDS 

We now address how boundary stacking can be used to en¬ 
hance the weak lensing signal of voids. Since it is a small 
effect, void lensing is difficult to measure (Melchior et al. 

2014) , although recently multiple detections of this signal 
have been reported (Clampitt & Jain 2015; Gruen et al. 

2015) . Increasing the signal to noise of this measurement, by 
either having a larger sample of voids and/or by improving 
how voids are stacked, would result in a powerful cosmolog¬ 
ical probe, especially for tests of modified gravity theories 
(Cai et al. 2015; Barreira et al. 2015). 

Within the thin lens and the Born approximation, the 
weak lensing signal is determined by the surface mass den¬ 
sity, 

X(Z)=P m J Sd,z)dz, ( 8 ) 

where £ is the position vector in the plane of the lens and z 
is the direction along the line-of-sight. We compute E(£) for 
three lines-of-sight that correspond to the simulation prin¬ 
cipal axes. For each line-of-sight we obtain E(£) on a 1280 2 
regular grid with grid spacing 0.39 /i -1 Mpc. We then pro¬ 
ceed to compute the lensing potential, 4/, via the relation 

V|tf(0*2^, (9) 

with the Laplacian operator restricted to the plane of the 
lens. The quantity, E c = c 2 Ds/I^tvGDlDls), is the criti¬ 
cal surface mass density for lensing, where Ds , 27l and Dls 
denote the angular diameter distance between the observer 
and the source, the observer and the lens, and the lens and 
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Figure 15. The stacked tangential shear, 7 t, of voids in three 
ranges in effective radius, ij e g. The top panel shows the spher¬ 
ically averaged result. The bottom panel shows the result when 
voids are stacked with respect to their boundary. The error bars 
show the l(j uncertainties due to object-to-object variation. 

the source. The exact value of E c , which depends on the 
characteristics of the lensing survey, is unimportant when 
comparing between the spherical and boundary stacking ap- 


proaches. 


For each point, we compute the convergence, 

re, and the 

shear, 7 = ( 71 , 72 ), as 


"(0 = 5[*11«) + * 22(01 = E (£)/S= 

( 10 ) 

71 (C) = \ [^n(C) - ^ 22 ( 0 ] 

( 11 ) 

72 (C) = ^ 12(0 = ^21 (C) J 

( 12 ) 


where the subscripts of T denote derivatives with respect to 
the two coordinate axes in the plane of the lens. 

For spherical stacking, the lensing signal is averaged as a 
function of the projected radial distance, r 2 D, from the void 
centre. This results in the convergence, re(r 2 D), which is a 
mean value inside a spherical shell of radius r 2 D ■ In the case 
of the shear, we are interested in the tangential component, 
7t, given by 

7 t = —71 cos(26>) + 72 sin(26 l ) , (13) 

where 9 is the angle between the first coordinate axis and 
the position of the point with respect to the void centre. Af¬ 
ter computing re(r 2 d) and 7t(r 2 D) for each void, we stack the 



rescaled projected radial distance r 9D / R eff 



Figure 16. Same as Fig. 15, but for the the stacked lensing 
convergence, re, of voids. 


voids according to their effective radius and across the three 
different lines-of-sight for which we computed E. Since the 
projected matter distribution is different along those orthog¬ 
onal lines-of-sight, averaging their lensing signal increases 
the signal-to-noise ratio. 

For boundary stacking, the procedure is slightly differ¬ 
ent, since we need to identify the boundary of the void in 
the lens plane. We do so by slicing the boundary of the void, 
which is a 2D surface, along the plane of the lens, with the 
slice centred at the point inside the void that is the farthest 
from the void boundary (this is the point corresponding to 
the minimum distance, V m in). Following this, we obtain a 
closed curve in the lens plane that corresponds to one par¬ 
ticular choice of the void boundary (see discussion below), 
which is then used to compute the distance in the plane of 
the lens, V 20 , of each surface element. Following this, for 
every void we compute the mean value of the convergence 
as a function of D 20 resulting in the quantity re(Z) 2 d). The 
tangential shear is computed using Eq. (13) but with 9 de¬ 
noting the angle between the first coordinate axis and the 
2D boundary distance vector, T> 2 d, at that point. Finally, 
we stack all voids of similar size and across the three lines- 
of-sight. 

We note that this is just one possible choice for stack¬ 
ing with respect to the void boundary, and may not be the 
optimal choice. For lensing studies, it is better to identify 
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2D voids in thin redshift slices, since this greatly enhances 
the lensing signal (Clampitt & Jain 2015). The boundary 
of these 2D voids is a ID curve in the plane of the sky. In 
such a case there is no ambiguity in choosing the ID void 
boundary in the plane of the lens. 

In Fig. 15 we show the void tangential shear obtained 
using the two stacking procedures. The spherically averaged 
7 t shows the characteristic dip of void lensing at r 2 D — Reff, 
which is nearly the same for the three void samples. This de¬ 
pression is more pronounced when using boundary stacking 
for which the signal is twice as large. Using boundary stack¬ 
ing increases the convergence, k, also by a factor of about 
two, as can be inferred from Fig. 16. This doubling of the 
lensing signal is the result of a better separation between the 
void border, where most of the mass is, and the void interior, 
which is mostly empty. This factor of two represents only a 
lower limit to the potential improvements resulting from the 
use of boundary stacking. Likely, the gain can be increased 
further by optimizing the selection of the void boundary in 
the plane of the sky. 


7 DISCUSSION AND CONCLUSIONS 

We have proposed a new method for characterising voids 
that has several advantages over the conventional spherical 
approach, as demonstrated by our analysis of galaxy voids 
in the Millennium cosmological simulation. This approach, 
which we call the boundary profile , is based on describing 
the structure of voids as a function of the distance from 
their boundary, which allows for a natural segregation of 
the inner, boundary and outer regions of each void. 

Voids are characterised by two defining features: they 
consist of large, fairly underdense volumes, with the evacu¬ 
ated matter found in a thin overdense region at the bound¬ 
ary, and they have very complex, non-spherical, shapes. The 
spherical averaging approach is inadequate for describing 
voids due to this very combination of features, as we ex¬ 
emplify for a simplified void model (Figs 1 and 2) and for 
realistic voids (Figs 4 and 5). This is a consequence of the 
fact that taking a spherical average over an intrinsically non- 
spherical object leads to a complex juxtaposition of the in¬ 
ner, border and outer regions of that object, with each re¬ 
gion having very different density. By contrast, the boundary 
profile method differentiates, by construction, between those 
regions. 

The boundary profile analysis revealed that the interior 
of voids is characterised by low densities that increase slowly 
towards the void boundary. This is followed by a steep rise of 
a density ridge at the void boundary, which decreases nearly 
as fast outside the void. The peak of the density ridge cor¬ 
responds to 1 + S ~ 4 while the interior of the void has 
1 + 5 ~ 0.2 — 0.4. We found a simple fitting function (Eq. 6) 
that describes fairly well the void density profiles and that 
can be parametrized in terms of a single quantity, the void 
effective radius (see Fig. 8). This parametrization provides a 
convenient way of describing the variation of density profiles 
with void size and allows for simple comparisons to theoret¬ 
ical models of void evolution, such as the spherical top-hat 
underdensity model (SvdW). 

The boundary density profile is self-similar, i.e. inde¬ 
pendent of void size, after rescaling the distance coordinate 


by the thickness of the void’s inner density ridge (see Fig. 
9). This suggests that the void interior knows about the void 
boundary or vice versa, and that the evolution of the two 
is coupled. This simple behaviour is reminiscent of the self¬ 
similar nature of dark matter haloes (Navarro et al. 1996, 
1997) whose origin, while not well understood, must reflect 
the scale-free nature of gravity. In contrast to haloes for 
which the characteristic scale is determined by the matter 
distribution in the innermost region, for voids the charac¬ 
teristic scale is determined by the matter distribution at the 
edge of the void. 

The boundary profile of the peculiar velocity reveals 
outflows from voids, which peak at a few Megaparsecs from 
the edge of the void, and an external infall region onto the 
void boundary. These outflows are preferentially directed 
along the radial direction, with the radial velocity being 
larger than the velocity component pointing towards the 
closest void edge (see Figs 12 and 13). The boundary profiles 
are especially suited for capturing the infall onto the void 
boundary, which is not seen for spherical profiles, and for 
determining if the voids are contracting or expanding (see 
Fig. 14). 

The boundary stacking method increases the weak lens¬ 
ing signal of voids by at least a factor of two when compared 
to the classical spherical stacking method. This gain can po¬ 
tentially be further increased by optimizing the selection of 
the void boundary on the plane of the sky (Cautun et al., 
in prep.). This gain in lensing signal boosts the utility of 
voids as cosmological probes, especially when applied to fu¬ 
ture large volume surveys like DESI (Levi et al. 2013), LSST 
(LSST Science Collaboration et al. 2009) and Euclid (Lau- 
reijs et al. 2011). 
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APPENDIX A: VELOCITY FIELD IN THE 
LINEAR APPROXIMATION 


In the linear approximation, the peculiar velocity at redshift, 
z = 0, is given by (Peebles 1980) 


Hf 

47rGp m S ’ 


(Al) 


where G is the gravitational constant and g is the gravita¬ 
tional field (for the remaining symbols see Eq. 7). The same 
relation holds for the velocity component, vii, along either r 
or T>, but with g replaced by g m . 

Applying Gauss’ theorem to the gravitational field, we 

have 




g ■ d S(x) = —4:nGM(< x ) , 
where x stands for either r or V and, 


M{< x) = 


nx 

J Xrr 


5{x')S{x')Ax , 


(A2) 


(A3) 


is the mass contrast enclosed by the surface, S{x), of con¬ 
stant x values. In the case of x = r, the lower integration 
bound, a: m m = 0, and the surface, S, corresponds to a spher¬ 
ical surface. For x = T), x m i n gives the distance, ©min, from 
the boundary of the farthest point inside the void (see Sec. 
3.2), while the surface, S, has an irregular shape. Fig. 1 
shows a cross section through an example void. Eq. (A2) 
can be rewritten as, 


_ AnGM(< x) 

511 " S[x) ’ 


(A4) 


where g« denotes the average value of g\\ over the surface 
^(a:). Inserting this last expression into Eq. (Al) results in 
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Figure Al. The area of a surface, S(T>), of constant void bound¬ 
ary distance, T>. It shows the mean area S(D), normalized by the 
area of a sphere of radius, T> + R ,\\, as a function of the normal¬ 
ized void boundary distance, T>/R e g. The three curves give the 
average value over all the voids in their respective R e g intervals. 


Eq. (7) used to compute vn lin (r) an d v\\ i in (X>). Note that 
since Eq. (7) is not linear in S(x), one needs to compute the 
linear theory predictions separately for each void, using their 
own density profile, and only in the final step to average over 
all the voids in the stack. 

To compute vm lin (0) one needs to know the function 
S(T>). This depends on the shape of the void boundary and, 
due to the large diversity of watershed void shapes, is dif¬ 
ferent for each void. Fig. Al shows the mean value of S(D), 
as measured for MS voids of various sizes. It shows that, 
when scaled appropriately, S(V) is approximately indepen¬ 
dent of the effective void radius. The scaled S(V) is maxi¬ 
mal for 0 = 0 since the void boundary is the most irregular 
T> = constant surface, as may be appreciated from Fig. 1. 
In the limit, V R e s, the surface S(T>) becomes a sphere 
and hence the scaled area shown in Fig. Al converges to 1. 
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